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We present a theoretical study of the physical properties of cationic lipid-DNA (CL-DNA) com- 
plexes - a promising synthetically based nonviral carrier of DNA for gene therapy. The study is 
based on a coarse-grained molecular model, which is used in Monte Carlo (MC) simulations of 
mesoscopically large systems over time scales long enough to address experimental reality. In the 
present work we focus on the statistical-mechanical behavior of lamellar complexes, which in MC 
simulations self-assemble spontaneously from a disordered random initial state. We measure the 
DNA-interaxial spacing, doNA, and the local cationic area charge density, gm, for a wide range of 
values of the parameter representing the fraction of cationic lipids. For weakly charged com- 
plexes (low values of <^c), we find that ^dna has a linear dependence on (^^^, which is in excellent 
agreement with x-ray diffraction experimental data. We also observe, in qualitative agreement with 
previous Poisson-Boltzmann calculations of the system, large fluctuations in the local area charge 
density with a pronounced minimum of cta/ halfway between adjacent DNA molecules. For highly- 
charged complexes (large <^c), we find moderate charge density fluctuations and observe deviations 
from linear dependence of dDNA on <j)c^ ■ This last result, together with other findings such as the 
decrease in the effective stretching modulus of the complex and the increased rate at which pores are 
formed in the complex membranes, are indicative of the gradual loss of mechanical stability of the 
complex which occurs when (j)c becomes large. We suggest that this may be the origin of the recently 
observed enhanced transfection efficiency of lamellar CL-DNA complexes at high charge densities, 
because the completion of the transfection process requires the disassembly of the complex and the 
release of the DNA into the cytoplasm. Some of the structural properties of the system are also 
predicted by a continuum free energy minimization. The analysis, which semi-quantitatively agrees 
with the computational results, shows that that mesoscale physical behavior of CL-DNA complexes 
is governed by an interplay between electrostatic, elastic, and mixing free energies. 



I. INTRODUCTION 



Somatic gene therapy holds great promise for future medical applications, for example, as new treatment for various 
inherited diseases and cancers [l|, Q • Viral vectors have been the most widely used systems for this purpose S 4 1 , 
but synthetic nonviral vectors are emerging as an attractive alternative because of their inherent advantages 2|- 
These advantages include ease and variable preparation, unlimited length of the transported DNA, and lack of specific 
immune response due to the absence of viral peptide and proteins [3 [1, Q . Complexes consisting of cationic lipids 
(CLs) and DNA comprise one of the most promising classes of nonviral vectors. They are already used widely 
for in vitro transfection of mammalian cells in research applications, and have even reached the stage of empirical 
clinical trials 1^ . Currently, their efficiency of gene transfer is considerably lower than that of viral vectors [ll|, ■ 
Substantial improvement of their efficiency is required before CL-DNA complexes become available for therapeutic 
purposes. 

CL-DNA complexes are formed spontaneously when DNA is mixed with cationic and natural lipids in an aqueous 
environment [3 US- Their formation is driven by the electrostatic attraction between negatively charged DNA 
and cationic lipid headgroups, and through the entropic gain associated with the concurrent release of the tightly 
bounded counterions from the CL and DNA [l^, [13, [H, flal. X -rav diffraction experiments have revealed that CL- 
DNA complexes exist in a variety of mesoscopic structures |T7l . ITsI . These structures include a multilamellar phase 
where DNA monolayers are intercalated between lipid bilayers (Lj^) [l^, and an inverted hexagonal phase with DNA 
encapsulated within cationic lipid monolayers tubes and arranged on a two-dimensional hexagonal lattice (Hjp) [l^ . 
In the more commonly observed L^^ phase, the DNA chains form a one-dimensional lattice, where the interaxial 
spacing doNA decreases with the charge density of the membrane. Isoelectric complexes, where the charges on the 
DNA exactly match those on the CL, are the most stable ones since they enable nearly complete counterion release 
p^ . For transfection, positively charged complexes are used, which can adhere to the negatively charged cell plasma 
membrane (T^. 

Despite all the promise of CL-DNA complexes as gene vectors, their transfection efficiency (TE; the ability to 
transfer DNA into cells followed by expression) remains substantially lower than that of viral vectors[ll|,[l2|- This has 
spurred an intense research activity aimed at enhancing TE i^lll, [l^, [H, [13, [HI [l^l ■ Recent tragic events associated 
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with the use of engineered adenovirus vectors have further stimulated the search for efficient synthetic DNA carriers 
[23l[2^. Recognizing that the structure of CL-DNA complexes may strongly influence their function and TE, much 
of the effort in theoretical and experimental studies has been devoted to understanding the mechanisms governing 
complex formation, structure, and phase behavior p^.[25|. The most widely used approach to describe the free energy 
of the complexes is the mean-field Poisson-Boltzmann equation which takes into account the electrostatic interactions 
of a charged continuum and the mixing entropy of the lipids and counterions [isl . [26l | . The bending energy of the lipid 
layer is introduced by means of the effective Helfrich energy p^. [27|. Thus, the preferred structure is determined by 
a delicate interplay between electrostatic interactions and bending elasticity, both depending on the molecular nature 
and composition of the lipids. The vast majority of neutral lipids, when complexed with CLs, lead to the L^^ phase. 
Addition of the neutral lipid DOPE (or other PE-based lipids) drives the spontaneous curvature negative, thereby 
inducing the transition to the Hjp phase . This transition is also promoted by addition of the cosurfactant hexanol 
which lowers the bilayer bending rigidity. 

The relationship between the physical properties and TE of CL-DNA complexes has been studied in a recent set of 
experiments utilizing a combination of several techniques (synchrotron X-ray diffraction for structure determination, 
laser scanning confocal microscopy to probe the interactions of complexes with cells, and luciferase reporter-gene 
expression assays to measure TE) [12]. The most notable result of these experiments is the identification of the 
membrane charge density, ctm, as a key universal parameter that governs TE of complexes The highest 

transfection rate has been observed at intermediate ctm, reaching values that are comparable to the high, ctm- 
independent, TE of DOPE-containing Hjp complexes [1^. That lamellar complexes compete the high TE of hexagonal 
complexes is of prime importance because, as mentioned above, most commonly used lipids prefer the Lj-" over the 
Hjp phase. Moreover, with newly synthesized multivalent lipids with headgroups whose charge is as large as Z = 5, 
it is possible to reach the optimal TE with a fewer number of CLs. This is a desirable feature which reduces the cost 
and, more importantly, the toxic effects of the CLs. It also means a smaller metabolic effort for the elimination of the 
lipids from the cell. 

Based on the above TE data, as well as X-ray diffraction and laser scanning confocal microscopy imaging, a model of 
cellular entry of L^^ complexes has been proposed which suggests that the process of transfection involves two stages 
p^ : (1) cellular uptake via endocytosis, and (2) escape of the complex from the endosome, presumably through 
the fusion of the complex with the endosomal membrane and release of the DNA into the cytoplasm. The adhesion 
of the complex to the cell is mediated by electrostatic attraction between the positively charged complex and the 
negatively charged cell's plasma membrane. The transfection process is limited by the rate of the second step, which 
increases exponentially with ctm- The independence of Hjp complexes TE on cta/ (in the low ctm regime) has been 
attributed to the mismatch between the positive curvature of the outermost lipid monolayer (which provides the 
complex with hydrophobic shielding) and the complex negative spontaneous curvature. This elastically frustrated 
state drives, independently of aM, the rapid fusion of the Hjp complex with the plasma or endosomal membrane. 
Lamellar complexes with very high ctm also exhibit reduced TE which should be attributed to the inability of the 
DNA to dissociate from the highly charged membranes (of the free complex) and to become available for expression 

Theoretical modeling of large molecular assemblies pose significant challenges due to the spatial complexity of such 
systems and by the range of temporal scales involved. In the case of CL-DNA complexes, the size of the complex 
may be as large as Ifim, while the basic unit cell of the L^^ complex is in the nanometer range (the DNA spacing 
is typically rfoNA ~ 20 — 70A, and the inter-layer spacing d ~ 65A). Since both short- (steric, hydrophobic), and 
long-range (electrostatic) interactions determine the physical and biomedical properties of CL-DNA complexes, it is 
essential that these systems will be studied at all possible levels of detail. Moreover, phase transitions of CL-DNA 
complexes as well as other topological changes (e.g., membrane fusion which occurs during transfection) involves the 
collective motion of many lipid molecules and, therefore, inherently take place on a variety of spatial and temporal 
scales. 

To address the multi-scale nature of CL-DNA complexes, a variety of models, differing in the length and the 
time scales of the phenomena of interest have been devised. At the microscopic molecular level, we have atomistic 
molecular dynamics (MD) simulations in which the lipids, DNA, and the embedding solvent are modeled explicitly in 
full (classical) atomic detail [s^l • These simulations provide valuable information regarding the molecular structure of 
the complexes, such as the role played by the neutral PC headgroups (more specifically, the N"*" end of the P~ — N+ 
dipole) in the screening of the electrostatic repulsion between the DNA chains. The length and time scales of atomistic 
simulations are limited by memory and CPU requirement to several nanometers and nanoseconds, which is far below 
the macroscopic regime encompassing the statistics and evolution of large molecular ensembles. At the macroscopic 
level, only the continuum behavior of existing CL-DNA structures can be addressed based on free energy functionals 
which are insensitive to the fine details of the lipids and DNA [l^, [2^, [26'| . Electrostatic screening effects between the 
DNA chains resulting from non-specific interactions between the lipids and DNA have been reported in these studies. 
This observation is complementary to the specific mechanisms observed in detailed atomistic computer models. 
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In this paper we present an "intermediate" molecular modeling approach that retains the most essential components 
of self-assembly and molecular statistics, but avoids the computational overhead of a full atomistic model. The 
model [sH extends an existing coarse-grained (CG) molecular bilayer model [H, [s^l by including both charged 
and neutral lipids, as well as charged DNA molecules. The inter- molecular potentials between different molecular 
species are designed to mimic the hydrophobic effect without the explicit presence of solvent. Thus, the approach 
carefully balances the need for molecular detail with computational practicality in a manner that allows for solvent- 
free simulations of complex self-assembly over long enough time scales to address experimental reality. In addition 
to showing spontaneous self-assembly of CL-DNA complexes, we also investigate the structural properties of lamellar 
complexes and measure a number of important quantities such as the dependence of the interaxial distance between 
DNA chains on the fraction of charged lipids, the polarization of the cationic charge distribution, local area density 
fluctuations, and the effective two-dimensional stretching elastic modulus, K*^. Some of these quantities are not easily 
accessible to theoretical continuum models and may require further approximations. One of the more interesting 
results is the decrease in K\ upon increasing the membrane charge density ctm, which reflects reduced mechanical 
stability and a higher probability of structural defects, such as membrane pores. The observations of such pores is 
consistent with the experimental findings of enhanced transfection efficiencies at high concentration of CLs, because 
pores must be formed to enable the escape of the DNA from the complex (see discussion above). It demonstrates 
the utility of CG modeling in addressing some key features of complex biological systems in general, and lipid-DNA 
assemblies in particular. 

We focus on isoelectric lamellar complexes in which the total charges carried by the CLs and DNA are equal. 
The mechanism of counterion release is most effective at this point, making the free energy of the complex minimal 
[l5| . The counterions concentrations inside the complex depend on their bulk concentrations. Here, we study the 
very low bulk concentration limit in which the counterions are (almost) completely depleted from the complex. This 
regime has been previously addressed in atomistic computer simulations of CL-DNA complexe s ISOll. Counterions 
effects have been dealt with in the framework of continuum Poisson-Boltzmann (PB) theory |25|. |26|. While this 
approach captures certain features of the charge distribution and the electric fields in the complex, it neglects several 
important factors such as the discrete nature of the ions and their finite sizes, which may govern their distribution in 
the gaps of (sub)nanometer dimensions that exist between the membranes and the DNA. Treating the water in these 
tiny voids as a bulk medium of effective dielectric constant e = 78 and neglecting dehydration effects are other gross 
approximations made by most current modeling techniques. 

We use a continuum model to analyze some of the simulation results. Our analytical study is based on the 
minimization of a phenomenological free energy functional with respect to the profile of the membrane and the 
cationic charge distribution. This free energy includes contributions from all the electrostatic interactions existing 
between the lipids and the (infinite array of) DNA molecules, as well as terms associated with the mixing entropy 
and (small length scale) protrusion modes [3J| of lipids. We show, both computationally and analytically, that the 
cationic charge distribution is polarized. The minimum of the charge density is obtained halfway between adjacent 
DNA molecules while the maximum is not reached right above the DNA, but is slightly shifted. The origin of this shift 
is the ability of the lipids not located right above the DNA to protrude and thus position their charged headgroups in 
regions of the complex where the electrostatic potential created by the DNA array is lower. Interestingly, we find that 
the charge density in the immediate vicinity of the DNA tends to match effective area charge density of the DNA rod. 
Charge matching between narrowly separated surfaces has been observed in previous studies of CL-DNA complexes 
[Tsl [23] , as well as in studies of other molecular assemblies [SSj . It has been attributed to the increased concentration 
of ions which are bound to remain in the confined volume between surfaces in order to neutralize the system. Charge 
matching is favorable because it enables the release of these strongly confined ions. Our study, which is based on an 
ion-free model, suggests that charge matching can be also driven by other mechanisms. 

The paper is organized as follows: In the next section we present our CG model of CL-DNA complexes and provide 
most of the details of the simulations (except for the details of a new Monte Carlo (MC) scheme that we use to sample 
the constant tension ensemble, which we introduce in Appendix [T]) . The results of our simulations are described in 
the two following sections dealing, respectively, with the dependence of the DNA spacing and the complex stability on 
the charge density, and the charge density fluctuations. The computational results are compared to the predictions 
of a analytical continuum model in Appendix [21 We close the paper with a brief discussion of the main results and 
an outline of some future prospects. 

II. COMPUTER MODEL 

The computer model of CL-DNA complex is based on a bilayer CG model presented elsewhere [H, HI] . The lipids are 
modeled as short trimer molecules consisting of one "hydrophilic" and two "hydrophobic" beads which are connected 
to each other by stiff linear springs The model does not include explicit solvent. Rather, a set of short-range 
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attractive intermolecular potentials is used, which effectively mimic hydrophobicity and allow self-assembly of bilayer 
from molecular disorder [3l|. Depending on the area density of the lipids, the bilayer is found in either a solid or a 
fluid phase [s^l, where the latter is characterized by an in-plane lipid diffusion and out-of-plane fluctuations whose 
spectrum is well depicted by Helfrich effective Hamiltonian 'SG^. The model of the CL-DNA complex is obtained 
by: (1) Choosing a fraction (pc of the lipids and placing a unit (point) charge -|-e at the centers of their hydrophilic 
bead. (2) Introducing DNA molecules, each of which is modeled as a rigid rod with a uniform axial charge density 
Adna = — e/1.7A and radius -Rdna = lOA. The diameter of the spherical particles constituting the lipids is set 
to cr ~ 6.3A (see definition of a in Ref. [s^P, which yields an area per lipid ajipid ~ IQK? for uncharged bilayers. 
Excluded volume interactions between rods (R) and spheres (S) are introduced via a truncated (at Tc = § + -Rdna) 

and shifted potential of the form: Ujis/kBT = 50 |[(cr/2 + -Rdna) /r]^^ — l|, where r is the distance between the 
center of the sphere and the axis of symmetry of the rod. The distance between nearest-neighbor rods is restricted to 
C^DNA > 2-Rdna- 

Modeling the DNA strands as infinite rods carrying uniform charge density A is consistent with the CG approach 
of the model, where only electrostatics, noise, and simple geometric features are retained. In this representation of 
DNA molecules we ignore the effects associated with their (1) flexibility and (2) the discrete nature of their charge 
distribution. The first approximation is justified in view of the fact that the DNA persistence length (^p ^ 500 A) is 
an order of magnitude larger than all the other relevant length scales in the problem. Curvature fluctuations involve 
free energy penalty of about 1 ksT per (,p of DNA length, which is negligible compared to the complex stabilization 
free energy of ~ 10^ — 10"^ ksT per persistence length [15|. The second approximation is supported by numerical 
studies revealing that the electrostatic potential around the DNA surface is not much different from that produced 
by the continuous charge density, except for a narrow regime in its immediate vicinity (37j . 

We study isoelectric complexes where the total charges of the DNA and the CLs neutralize each other, with no 
added counterions. Simulations of the quasi two-dimensional (2D) complex are conducted in a rectangular system of 
size Lx X Ly X Lz, with full periodic boundaries along the x and y directions, and periodicity with respect to only lipid 
mobility and short-range interactions in the z direction. The simulations were performed at room temperature and with 
a bulk water uniform dielectric constant e = 78. The rods are arranged in a ID array, parallel to the y axis and with 
equal spacing along the x-direction. Long range electrostatic interactions between the charged spheres, infinite rods, 
and their periodic images were accounted for using the Lekner summation method ;38]. The electrostatic potential 
energy, per simulation cell, between a CL whose charged headgroup is located at f+ Ar = {x + Ax, y + Ay, z + Az) 
and another CL and its replicas located at {x + niLx, y + nLy, z), where m, n are integers, is given by the exponentially 
convergent summation of, e.g., the form: 



.A^^ e2 f 4 ^ / Ax \ 
Vss (Ar) = -|^E-«(2.-nj 



E ^0 



fc=-c 



27rn 
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, , „ Az\ / Ay 

cosh I 2tt I — cos I 27r^— 



(1) 



where Kq is the modified Bessel function of order 0. The self-energy V^g that arises from a charged sphere with its 
own periodic images is found by evaluating the expression 



(2) 



and is given in Ref. [38| . The sphere-rod electrostatic energy per simulation cell is the combined logarithmic interac- 
tions between a point charge and a ID array of line charges. It is given by f39l| 
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(3) 



The rod-rod electrostatic self-energy is [3£ 



(4) 
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Simulations of electrostatics in water-free models are usually conducted with a uniform dielectric constant due to 
the complexity of including mirror charges in a disordered molecular ensemble. In PB theories [Tsl . [35| a different 
type of boundary conditions (BCs) is usually assumed, namely that the dielectric constant vanishes in the interior of 
the DNA and the lipid membrane. This latter approximation is justified by the fact that the (bulk) water dielectric 
constant, e — 78, is much larger than all the other relevant dielectric constants. The electric field lines prefer to stay in 
high dielectric media and the limit e = corresponds to systems in which the electric field is entirely contained within 
the aqueous part. Fortunately, the exclusion of electric fields from the DNA and the membranes is not solely related 
to their low dielectric constants. In our simulations, it is mainly a matter of the geometry and the charge distribution 
in the system. Therefore, the electrostatic forces, through which the charges in water interact with each other, are 
expected to be insensitive to the dielectric constants of the DNA and the membranes. Computational studies of 
electrostatics near similar simple geometric interfaces indicate that the net effect of mirror images is, indeed, minor. 
This observation also serves as a justification for our choice of BCs. Rather than using full periodic BCs in all three 
directions, we study an infinite "slab" consisting of a single array of charged DNA molecules and two bilayers (see 
Fig. H]), in which only the inner monolayers (facing the DNA rods) are charged while the outer monolayers consist of 
neutral lipids only. The interaction between different slabs (whose overall net charge is neutral) is well screened and 
is significantly weaker than the Coulomb interactions between the charged components within each slab. 

In a preceding publication [3l| . we have shown that complexes such as the one appearing in Fig. [T] are formed 
spontaneously in simulations starting from a disordered initial state where the lipids are randomly distributed within 
the simulation cell. This demonstrates that the complex represents a stable equilibrium phase of the system. In 
this work we focus on structural properties of CL-DNA complexes and use pre-assembled complexes for this purpose. 
The simulations were performed at constant surface tension 7 = by employing a new sampling scheme to generate 
area-changing trial moves. The sampling scheme, which is different from the commonly used method of sampling the 
(iV, 7, r) ensemble [1^ , is described in detail in Appendix [T] The rest of the details of the simulations appear in 
Ref. [3l|. 



III. DNA SPACING AND MECHANICAL STABILITY 



Measuring the DNA interaxial spacing, doNA, serves as a critical test to our model's ability to mimic the mesoscale 
behavior of CL-DNA complexes, because doNA can be measured in X-ray diffraction experiments. The experimental 
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data of Safinya et al. [13|, |18| shows that for isoelectric complexes, the dependence of ^dna on the membrane charge 
density ctm — 2e(/)c/aiipid is governed by the relationship 

dDNA = = ( J (5) 



which results from simple mass conservation in the lamellar geometry. 

The computational results for the average spacing between adjacent DNA rods, doNA, are plotted in Fig. [2] as a 
function of the inverse of the fraction of charged lipids, l/<^c- The solid line is a fit to Eq. ([5]) with anpid = 69A^ which 
is the area per lipid in uncharged membranes. The deviation from linear behavior at high charge densities arises from 
the increase in aiipid with (^^ (see Fig. [3]). The numerical data is in excellent agreement with the experimental results 
reported in Ref. [13] ■ Specifically, the experimental ^dna vs. 0c data [see Fig. 4 (B) in Ref. show agreement 
with Eq. ([5]) at low charge densities (with a value of aupid which is slightly different than the one defined here) and a 
similar deviation trend at high charge densities (note that our Eq. (O and the comparable one in Ref. ^Tl| express the 
same relationship in different forms. See discussion in [4l[.) The assumption underlying Eq. ^ is that the effective 
interactions between the DNA are repulsive and balanced by the elastic membrane forces. Linear elastic stress acting 
on a membrane is related to aupid and its equilibrium value aj'jpjj by t = Ka (^aiipid — '^.^pj^^ /'^lipid' where Ka is the 

2D stretching modulus, which for lipid bilayers is typically in the range Ka ^ 10^ ergs/cm^. At high charge densities, 
the electrostatic stress is sufficiently large to eliminate the membrane thermal undulations and increase aupid (4^ . In 

the present study, we find for complexes with ^ 0.85 that the strain e = ^^aupid — aiipid) /"^npid ^ ^■'^1 which is 
somewhat larger than the typical strain lipid bilayers withstand before rupture (e 0.02 — 0.05, [il|). The discrepancy 
between experimentally observed rupture strains and the strain found in our model is not surprising given the coarse 
grained lipid-model of simple three point objects. It may also be partially attributed to the system size dependence 
of the rupture strain p4| . Membranes with higher have indeed been found to be susceptible to pore formation, 
as illustrated by the configuration in Fig. [5] of a complex with <^c ~ 0.9. The loss of mechanical stability is also 
evident from the rapid decrease in the effective stretching modulus of the complex K\ for 0c ^ 0.7 (Fig. 2]), which 

has been extracted from the mean square of fluctuations in anpid: K\ — ^bT 



N (aiipid - aHpid) 



The 
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larger area fluctuations at high <j)c increase the probabihty of pore opening which in turn may lead to disassociation 
of the complex. 

It is interesting to compare the inter-DNA electrostatic interactions in the complex with the same interactions 
in bulk in the presence of monovalent counterions ^45*1. In both cases, the bare electrostatic repulsion is screened 
by the distribution of microions and CLs around the DNA. In the complex, the distribution of the cationic charge 
is limited by a geometric constraint, namely the residence of the CLs on the membrane surface which maintains a 
finite separation from the inter-DNA plane. Therefore, we expect that screening to be less efficient than in bulk 
solutions. It is also reasonable to expect that, compared to 3D systems, the confinement to a surface increases the 
repulsive electrostatic interactions between the cationic charges which further increases the magnitude of the effective 
inter-DNA repulsive force. 

Since our model does not include water explicitly, the hydration forces which may be significant at small cJdna 
are missing from the picture [1^. This can be corrected either by introducing an additional short range DNA-DNA 
potential that explicitly account for the hydration free energy or assuming that i?DNA represents the hydrated rather 
than the bare DNA radius. Since the hydration forces decay on DNA-DNA surface separations of more then ~ 1 nm, 
their introduction into the model will only modify the results at small c^dna (large (pc)- Specifically, this will increase 
the effective DNA-DNA repulsion and, therefore, will strengthen the trend observed in Fig. [2] that at high charge 
densities Eq. ([5]) underestimates doNA- It will also shift the limit of mechanical stability to slightly lower charge 
densities. Another feature missing in our model is the contribution of DNA bending fluctuations to the inter-DNA 
interaction. The effect of these fluctuations in bulk is to increase the decay length of both hydration and electrostatic 



forces [47[. One may expect a similar contribution to inter-DNA interactions in CL-DNA complexes, although we 



are unaware of any systematic study of this effect in 2D. X-ray studies find weak positional disorder in CL-DNA 
complexes. The typical correlated domain size of the ID lattice of DNA extends to nearly 10 unit cells [l^, which is 
twice as large than the size of the complex in our simulations. 

Our computational results explain well the recently observed enhanced TE of lamellar CL-DNA complexes at high 
charge densities [H, HH] . The limiting stage of the transfection process is the escape of the complex from the endosome 
in which it is initially trapped after entering the cell. Escape from the endosome occurs through activated fusion of the 
complex and endosomal membranes, during which both must be perforated. Having a complex with poor mechanical 
stability is an advantage at this stage, since such a complex will tend to open pores more easily, and through these 
pores the DNA may be released to the cytoplasm. We suggest that the loss of mechanical stability results from the 
cationic charge of the lipids and the pressure which it exerts on the complex membrane. At high charge densities this 
pressure exceeds the rupture tension of the membrane and, thus, leads to mechanical failure of the complex. 



IV. CHARGE DENSITY MODULATIONS 



The quantity defined as <f>c represents the mean number fraction of charge lipids. However, the lateral distribution 
of cationic charge on the membranes need not be uniform. One may expect the CLs to accumulate above and below 
the negatively charged DNA rods. This tendency to minimize the electrostatic energy of the CLs-DNA interactions 
is competed by the thermally induced mixing entropy and the repulsive electrostatic interactions between the CLs, 
which favor homogeneous composition of the cationic and neutral lipids. Furthermore, charge density modulations 
may be coupled to membrane undulations p?! . |48| , and both can contribute to lowering the free energy of the complex. 
As discussed above, the stability of the complex is directly related to its TE and, therefore, it is important to study 
the effect of these "degrees of freedom" of the lipids. 

The dependence of <f>c on x, the position within a unit cell of the complex (i.e., the interval between adjacent 
DNA rods) is depicted in Fig. [51 where x = Q and x — c?dna correspond to lipids located right above or below the 
DNA (see inset of Fig. (5]). The curves, from bottom to top, are for the following value of the mean number fraction: 
(j)^ = 150/500 = 0.3; 150/375 = 0.4; 150/300 = 0.5; 150/250 = 0.6; 150/220 - 0.68; 150/195 ~ 0.77; 150/185 ~ 0.81. 
As expected, we find that for all values of the minimum of 0c(a^) is achieved for x = c?dna/2, i.e., in the middle of 
the unit cell. The minimum is more pronounced for low values of 0c, in which case the maximum of (j^ci^) is at a; = 
and X = rfoNA- At the higher values of (j>c, the maximum shifts from the edge of the unit cell towards the center and, 
in general, the fiuctuations in 4>c{x) become quite small. 

The shift in the maximum of (pdx) from the immediate vicinity of the DNA towards the center of the cell has been 
previously reported in theoretical studies of the system based on PB theory [Tsl . l27j . It has been attributed to the 
tendency of the system to match the areal charge density of the membrane with the effective areal charge density of 
the DNA, (Tdna = — ADNA/(27ri?DNA) ~ 9.4 • 10~'^e/A^. This involves attraction of CLs towards the DNA at low 
values of 0c and significant charge modulation over the relatively large distance between the DNA rods. On the other 
hand, when 0c is large and doNA is small, the charge density fiuctuations are weak and CLs must be depleted from 
above/below the DNA to match the local charge density of the DNA. 
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The tendency to match the local charge densities of the membranes with ctdna is seen more clearly in Fig. [T] Here, 
the charge density 17^(2;) rather than (j)cix) is plotted as a function of x. The dashed horizontal line corresponds 
to the effective charge density of the DNA ctdna ~ 9.4 ■ The graphs show that for all values of (j)c, 

the local charge density at the edge of the unit cell remains within 15% of ctdna- More interestingly, the graphs 
show that in contrast to (jjcix), the maximum of (7{x) is always shifted from the DNA towards the center of the 
unit cell. This observation deserves some special comment; In most analytical studies of membranes, area density 
fluctuations are neglected (or, rather, it is assumed that the local area per lipid aiipid(a;) is constant) and, therefore, 
4>c{x) and aM{x) = 2e(f)cix) / anpidix) are proportional to each other. Our computational results indicate that area 
density fluctuations (see Fig. [8]) may be quite important and serve as an additional degree of freedom that further 
reduces the free energy of the system. In Fig. [51 po = ('^Upid)^^ — 1/69A~^ denotes the area density of uncharged 
membranes. For large (j)c, P < Po for all values of x, reflecting the fact the mean area per lipid increases at high 
charge densities (see Fig. [3] and discussion above) . More noticeable are the area density fluctuations within the unit 
cell which can be observed for all values of (f>c- The location of the maximum area density coincides with that of the 
maximum charge density and, therefore, can be attributed to the accumulation of charged lipids. Had the area density 
been constant, this would mean depletion of the neutral lipids from the same region in the unit cell. Area density 
fluctuations represent an additional degree of freedom of the system, which permit a more uniform distribution of the 
neutral lipids, and, thus, pays off in terms of lower mixing entropy. The magnitude of the area density fluctuations is 
roughly given by y/ {p/po — 1)^ ~ (^bT/ K'^^ay^^-^^y^'^ , which for typical values of the parameters {K'^ 250 ergs/cm^, 

^Upid ^ 70A^ - see Figs. |3] and HI yields \/ [p/po — 1)^ ^ 0.15, in reasonable agreement with our results in Fig. [S] 

As mentioned above, local charge density matching is the origin of the CLs tendency to migrate towards the middle 
of the unit cell. Solutions of the PB equation [l^, [s^ show that the concentration of counterions, which will be bound 
in the narrow water gap that exists between the DNA and the membrane, increase with the charge density mismatch. 
An accumulation of counterions in such a small volume is energetically unfavorable and will lead to a very large 
osmotic pressure. Two comments should be made regarding the solution of the PB equation: First, the screening of 
the electrostatic interactions by the highly confined counterions is probably overestimated by continuum PB theory 
because the slithering of counterions into the small gaps will be hindered both dynamically (excluded volume) and 
thermodynamically (dehydration). Second, our simulations of isoelectric complexes with no counterions apply to 
the no-screening limit where this effect is not expected to occur. We therefore conclude that local charge density 
matching may be also driven by other factors. Indeed, in the absence of screening one must consider the interactions 
of the CLs with the periodic array of line charges rather than with the closest DNA rod. The Coulomb energy due 
to the interaction of a charge +e with an infinite periodic array of rods of charge density per unit length A < is 
given by Eq. ([3]), with L^. = c?dna- For a charge residing on a perfectly flat surface located a distance Az = D = 
-RDNA+cr/2 ~ I3A above the mid-plane of the DNA array, Eq. ^ simplifies to Vrs (Ar) ~ (eA/e) i?cos (27rAx/c?DNA), 
with B = 2 exp (— 27r_D/c?DNA)- This expression is valid as long a,s B < 1, which is indeed the case for the above value 
of D and the range of values of the DNA interaxial spacing duNA ~ 25A— 50A considered in this work. For a nearly 
flat surface with Az = D — h {0 < h <^ D), the sphere- rod electrostatic energy reads 

, . . eA r 27r/i / 27rA2;\l , . 

Vrs{A^-— -, -fScos . (6) 

e L"DNA V"DNA/_ 

The first term in this equation reflects the long-range nature of unscreened electrostatic interactions which makes the 
potential of the infinite array of line charges look similar to the potential of a uniformly charged surface with areal 
charge density um — A/doNA- The second term is due to the periodicity of the system and represents the tendency 
of cationic lipids to favor the proximity of the anionic DNA rods {Ax = 0). The attraction of the CLs to the DNA 
rods, located at the edge of the unit cell, will be offset by their attraction towards the mid-plane of the DNA array 
[first terms in Eq. ([6])]. This attraction draws the CLs towards the center of the unit cell because, right above the 
DNA rod, the vertical separation between the surface and the DNA array is restricted to Az = D {h — 0) by excluded 
volume interactions. Away from the DNA and close the center of the unit cell, the elastic deformation of the surface 
permits Az < D {h > 0) which is energetically favorable. The highest cationic charge density will be obtained at the 
minimum value of the sphere- rod electrostatic energy Vrs- A detailed calculation based on a continuum expression 
for the free energy which includes contributions of the electrostatic and elastic energies and of mixing entropy, is 
presented in Appendix [2l We show that only for infinitely rigid surfaces the maximum of charge density is observed 
at the edge of the unit cell. In all other cases, the membrane tends to deform towards the mid-plane of the DNA 
array, which leads to shifting of the maximum of <tm{x) towards the center of the cell. This is in agreement with our 
results in Fig. [71 although the deformation of the membrane is unnoticeable in snapshots of the system (e.g., Fig. [T]). 
Our calculation shows that the typical amplitude of the deformation is extremely small, of the order of 1 — 2A. This 
estimate is model dependent, but in agreement with previous studies of the system [13,14^, and explains the apparent 
flatness of the membrane observed in our simulations. 
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V. DISCUSSION 

We have presented a molecular simulation method which captures the self-assembly of cationic liposomes complexed 
with DNA - a promising synthetically based nonviral carrier of DNA for gene therapy. The method is an intermediate 
modeling approach between atomistic computer simulations and continuum phenomenological theories. Like the 
former, it utilizes a molecular description of the system; but similarly to the latter, it employs a coarse-grained (CG) 
representation of the intra-molecular atomic details. The reduced number of degrees of freedom, as well as the fact 
that the model does not require explicit representation of the embedding solvent, lead to a significant improvement 
in computational efficiency. Thus, the approach carefully balances the need for molecular detail with computational 
practicality in a manner that allows for solvent-free simulations of complex self-assembly over long enough time scales 
to address experimental reality. 

In addition to showing spontaneous self-assembly of cationic lipid-DNA complexes, the broad utility of the model is 
illustrated by demonstrating excellent agreement with X-ray diffraction experimental data for the dependence of the 
interaxial distance between DNA chains, dDNA7 on the fraction of charged lipids (t)c- Specifically, we find that ^dna 
is inversely proportional to (pc - sl relationship which can be also derived by a simple packing argument where the 
DNA rods form a space-filling ID lattice. This result is indicative of a repulsive long-range inter-DNA interaction. 
The predominant contribution to this interaction is due to non-specific electrostatic repulsion between the negatively 
charged DNA rods, which is only partially screened by the cationic charge on the membranes. We note that the 
magnitude of the repulsive interaction plays no role in the packing argument. Therefore, a linear relationship between 
doNA and cj)'^^ is predicted by both our CG simulations and the PB theory, despite of the fact that in the PB treatment 
the screening of inter-DNA repulsion is also due to the counterions presented in the complex. 

Certain features of CL-DNA complexes, for instance the process of self-assembly and structural defects (Fig. [5]), 
can be addressed more effectively through CG simulations rather than by continuum theories of existing structures. 
This point is nicely demonstrated by our simulations of highly-charged complexes. Upon increasing the fraction of 
the CLs, we find that: (1) the area per lipid increases, (2) the effective stretching modulus of the complex decreases, 
and (3) the rate of pore formation increases which eventually leads to the disintegration of the complex. All together, 
these results indicate that the higher the charge density of the membranes the lower the mechanical stability of the 
system. This is a key observation which may explain the recently observed enhanced transfection efficiency (TE) 
of lamellar CL-DNA complexes at high charge densities. Transfection is viewed as a two-stage process: (1) cellular 
uptake via endocytosis, and (2) escape of the complex from the endosome, presumably through fusion of the lipids 
with the endosomal membrane and release of the DNA into the cytoplasm. TE of lamellar complexes is limited by the 
rate of the second stage and, hence, increases with the decrease of mechanical stability, i.e., with increase of charge 
density. 

Given the consistency of agreement between our CG molecular approach and observed experimental features, we 
suggest that the presented model is an appropriate and promising tool for investigating the statistics and dynamics 
of lipid-DNA complexes on spatial and temporal scales relevant for biological and biomedical applications. In future 
work, we plan to develop models which would mimic CL-DNA complexes with improved gene delivery performance, 
such as complexes containing multivalent lipids and lipids attached to short polymer chains. A special effort will 
be made to develop a model for the inverted hexagonal (Hjp) structure, and to examine the mechanical behavior 
of this phase, which appears to be experimentally quite distinct from the behavior of the lamellar phase. We will 
also investigate the effect of counterions which must be presented in the positively charged complexes that adhere 
to the negatively charged cell membrane (at the initial stage of the transfection process). The model may be also 
extended to include some features of the DNA helical structure. These more advanced models may lead to a better 
understanding of the principles governing the statistical-mechanical behavior of CL-DNA complexes, which is crucial 
for systematic and successful design of efficient synthetic vectors for gene therapy. 

APPENDIX 1: SIMULATIONS AT CONSTANT SURFACE TENSION 

Bilayer Membranes may be considered as narrow interfaces consisting of lipids and the hydration layers that 
separate two aqueous bulk phases. Simulations of liquid/liquid interfaces can be performed in a variety of statistical 
ensembles. Assuming that the temperature T and number of particles A'^ are fixed, one may use the (TV, T, V, Ap) 
ensemble, in which the total volume of the system, V, and the projected area of the interface, Ap, are held constant. 
Alternatively, the normal and transverse components of the pressure tensor, P„ and Pt may be fixed, letting V 
and Ap fluctuate. Another common choice is the constant surface tension ensemble {N,T,V,"f), which mimics the 
experimental conditions more closely than the {N, T, V, Ap) ensemble does. The {N, T, V, 7) is of particular importance 
for simulations of membranes, which can exhibit large undulations at vanishing surface tension. Accessing the 7 — 
regime is crucial to modeling such systems. 
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There has been an ongoing theoretical debate concerning the statistical thermodynamic definition of 7. In analytical 
studies the surface tension is usually regarded as the thermodynamic variable conjugate to the total interface area 
A, which is the sum of the projected area Ap, and the area stored in the thermal undulations A A. Because of the 
relatively high value of their stretching modulus, bilayer membranes are often assumed to have a fixed total area, and 
7 is used as a Lagrange multiplier fixing the value of A. However, in computer simulations it is difficult to sample 
an ensemble where A is constant since AA cannot be easily controlled and, moreover, its value is not well defined 
(in contrast to continuum models). It is, therefore, more common in computer simulations that the surface tension is 
treated as conjugate to Ap, which is the cross-sectional area of the simulation cell. With this convention, one readily 
derives the following relationship between 7 and the pressure tensor [ioj 



7 = (L„ X (P„ - Pt)) , 



(11) 



where L„ is the size of the system in the direction perpendicular to the interface, and () denotes thermal average. 
This quantity coincides with yet another quantity commonly referred to as the "surface tension" , namely the 
coefficient in the expression describing the dependence of the mean thermal fluctuations on the wave vector (the 
"spectral intensity") [ssj : 



knT 



aiipid [jq^ + Kq^ + Oiq^)Y 



(12) 



To simulate the (N, T, V, 7) ensemble, one needs to sample configurations, in which the total volume of the system 
is conserved while the area is allowed to fluctuate. The common method to generate such an ensemble is to consider 
a rectangular simulation box of volume V — x Ly x and projected area Ap — x Ly and, occasionally, rescale 
the dimensions of the box and the molecular coordinates {r *} in the following manner [40|: 



Ly > Ly + SLy ; Vy 



Lx + 5Lx 

L'T. 



Ly + SLy 



L. 



V 



{Lx + SLx){Ly + 5Ly) 



LxLy 



{Lx + SLx){Ly + SLy) 



(13) 



In scaled coordinates, {/ * = (r^/ia;, /Lj,, /L2)}, the partition fmiction is given by 



Z = 



dLx dLy dLz 6 L^ 



(l, - ■ f nili< dii dii v^e-P^^^ g-/3a({r'},L.,L„L.)^ (^4) 

V LxLy) Jq 



where 6 is the Dirac delta function, U is the energy of the configuration and /3 — l/ksT. Since the volume V and 
number of particles N are both fixed, the acceptance criterion is given by 



-Pacc(o n) = min 1, e 



(15) 



where 6U and 6Ap denote, respectively, the difference in the energy and area between the new (n) and old (o) 
configurations. 

Eq. (jl3p defines a one-to-one, locally volume preserving, transformation of the molecular coordinates between 
rectangular simulation boxes of slightly different dimensions. In principle, however, only the total volume of the 
simulation cell must be conserved and, therefore, other one-to-one transformation may be proposed. One such 
transformation, which is particularly suitable for solvent-free interfacial systems, is the following: 




(16) 



(L^+SL^){Ly+SLy) 

mod(r^ + (5r 1,^2(0) )" 



Sri 



i = 1 
1^1 



where Lz(0) is the initial (at t = 0) height of the simulation box, and "mod" is the modulus operator. Unlike 
transformation only the coordinates of one particle [labeled "1" in Eg. p?)) ] are rescaled proportionally to the 
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size of the box in (fT7|) . Therefore, a-priory, only this particle is guaranteed to remain within the rescaled simulation 
box, while all the others may be displaced beyond the boundaries of the system in the z direction (0 < z < L^). 
Detailed balance requires that in the case of such an event, the move attempt will be rejected, i.e., the energy of the 
new configuration is defined U — oo. Let us consider an interface located close to the center of the box (z ~ Lz/2), 
separating two dense bulk liquid phases. Let us also assume that the particle with index "1" resides on the interface. 
For such a solvent-containing system, rescaling the dimensions of the simulation box and the molecular coordinates 
according to transformation P7)l will usually fail. Decreasing Lz [Lz < ^^(O)] will lead to ejection of some solvent 
particles from the system, while increasing Lz [Lz > Lz{0)] will lead to the formation of an empty stripe and changing 
the bulk densities which is energetically very costly. In solvent-free models, on the other hand, the bulk phases 
are "empty" and such a problem will arise only if the interface is located close to one of the z boundaries. For 
infinitely long runs, the fraction of time that the interface spends near the boundaries scales like w/Lz (where w 
is the physical width of the interface), which can be made arbitrarily small by increasing Lz (i.e., by changing the 
"volumes" of the empty bulk phases). In practice, Lz need not be very large since the diffusion of the interface is 
so vanishingly slow that it will never reach one of the boundaries within conceivable simulation time. With function 
(|17p defining the transformation between simulation boxes of different shapes, and in terms of the scaled coordinates 
{I ^ — {rl/Lx,ry/Ly,rl/Lz);l = (r^/Lx, r* /Ly, r* /Lz(0))}, the partition function in (fT4|) is rewritten: 

f°° f V 

Z = I dLx dLy dLz S [ Lz 

Jo \ LxLy 

°\tidlldlldllV[Lz{Q)Aj,f-^ ^-PiA, ^-pu{{i^}x^,L,x.,LAo)) ^ ^^7) 



from which we readily conclude that the acceptance criterion should be 



:(0 



1, 1 



5A. - 



Ap 



P \ ^-l3{5U+'ySAp) 



(18) 



This criterion resembles the acceptance criterion for simulations of the 2D isobaric-isothermal ensemble, accept for 
the exponent of the term (1 -I- 6Ap/Ap) being — 1 rather than N. This could be interpreted in the following 
way: Simulating a thin interface of width w at constant 7 is similar to simulating a 2D system at constant pressure. 
The center of mass of the interface can be found with equal probability at any position along the z direction of the 
simulation box. However, because of the total volume conservation, the height of the (mostly empty) simulation box 
scales with the inverse of the area Ap. The contribution of this degree of freedom is expressed by the additional factor 
il + SAp/Ap)-\ 

We have performed simulations of identical systems (both neutral membranes and CL-DNA complexes) using 
sampling methods (|13[) and (jl7p . and measured the probability distributions of the energy and the area of the system. 
Both methods produced identical distribution functions, but the time it took to obtain reasonably accurate results was 
considerably shorter with the new sampling scheme [Eq. ^T7\ ] (about 2-10^ MC time units) than with the conventional 
one [Eq. (roughly lO"^ MC time units). The superior effectiveness of the new scheme should be attributed to 

larger area changes per reshaping attempt, SAp, that the new scheme permits, which were typically half an order 
of magnitude larger than in the old scheme. We believe that the origin of this is the fact that the membranes in 
our simulations are "softer" with respect to area changes than they are with respect to the variations of their width. 
In the conventional sampling scheme, area and width fluctuations are coupled by local volume conservation, which 
greatly reduces the magnitude of acceptable reshaping moves. In the new scheme this coupling is removed, enabling 
area variations which do not simultaneously squeeze the layers against each other or pull them apart. 



APPENDIX 2: ANALYTICAL MODELING 



In order to estimate the deformation of the membranes and the charge density fluctuations within a unit cell of the 
complex, we consider the system shown schematically in Fig. (S) This system consists of three charge distributions: 
an infinite array of equally spaced rods with density per unit length A < 0, and two surfaces with mean charge 
density ctm > per unit area (representing the monolayers facing the DNA array on each side). The vertical distance 
Az{x) of the surfaces from the mid-plane of the DNA is Az — D = -Rdna + (7/2 ^ ISA at the edge of the unit cell 
[x = 0, ^dna), and may be smaller [Az ^ D — h < D ; h{x) > 0) for < x < ^dna due to the electrostatically 
induced deformation of the surfaces. We denote the charge density fluctuations by 6(j{x), and consider the limit of 
small fluctuations \Sa\ <^ om, as well as the limit of small membrane deformation |/i|/(iDNA ^1. As a reference state 
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for energy calculations, we take a complex with perfectly flat surfaces and no charge density fluctuations. Since the 
complex is isoelectric, the mean charge density of the surfaces is related to the DNA linear charge density by 



2a M 



X 



(^DNA ' 



while 



^0 



icTM {Vh{x)f + Sa (x) 



dx 



da (x) dx = 0. 



(21) 



(22) 



For a fixed value of dDNA, the free energy of a unit cell of the complex per unit length in the y direction (parallel 
to the DNA rods) consists of the following contributions: 

1. The Coulomb energy of the interaction between the DNA rods and the charged surfaces which is given by [see 
Eq. ®] 



[aM + 5a{x)] 



A 



2'Kh (x) 



B cos 



/ 2t:x \ 

\ C^DNA / 



l + \{^h{x)f 



dx. 



(23) 



The prefactor 2 in this expression is due to the two surfaces interacting with the DNA rods. Treating B in Eq. 
as a small parameter [B = 2 exp (— 27r_D/(iDNA) ^ 0.12, for D ^ IsA and cZdna ~ 30A, which is of the same order of 
magnitude as (da /a) and (^/rfoNA)]; we see that the leading term in the above expression is 



F 1 = 2 



aM- 



X 2Trh (x) _ _ 4o-^,/rfDNA 



e rfoNA 



2'Kh (x) 

C^DNA 



dx. 



(24) 



The superscript in (as well as in the other terms of the free energy appearing below) denotes the order of the term 
in the small parameters of the expansion. The next (second order) term is given by 



Sa{x) 



A/ "DNA 



27r/i (x) 

doNA 

I A 

Sa{x) 



B cos 



/ 27rx \ 

V rfoNA / 



2Trh (x) 



^DNA 



B cos 



dx 
/ 2ttx \ 



dx. 



(25) 



2. The Coulomb energy of the interaction between the two charged surfaces, which in leading order is given by 



F?^ 



5a{x) — ; — ^^-^dx. 



IDNA 



(26) 



3. The Coulomb energy of the interactions between lipids residing on the same surface (the surfaces self electrostatic 
energy) [see Eq. ^ with A [a^f + Sa{x')]dx' and Az/L^ —h/doNA 0, and recall Eq. ([^ ] 



F' 



dx ■ 



O-M^DNAg f h{x) 



IDNA 



2 - 2 cos 27r 



doNA 



(27) 



where G (/i/c^dna) is a dimensionless function. 

4. The mixing entropy of the charged and neutral lipids in each surface which, ignoring the variations in the area 
per lipid (see Fig. [5]), is given by 



knT 



Sa^ (x) 



"lipid Jo CTAf (e/fliipid — (TAf) 



dx. 



(28) 



where e/aiipid is the maximum possible charge density (obtained when all the lipids are charged) and, obviously, 
ctm < e/oiipid. 

5. The elastic energy of the deformation of the surfaces. In previous theoretical studies of CL-DNA complexes 
[27l . IZsj , this energy has been associated with the bending of the surfaces and has been expressed in terms of Helfrich 
effective Hamiltonian [sHi- However, Helfrich effective Hamiltonian captures the elasticity of surfaces only on length 
scale which are typically larger than the length of the unit cell ^dna ~ 25 — 50A. At smaller scales, the elastic behavior 
of membranes is dominated by individual or collective lipid protrusions. Protrusion modes tend to increase the local 
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surface area and the restoring force acting against these protrusions can, therefore, be characterized by an effective 
local surface tension 7^. The corresponding elastic energy (per unit length in the y direction) is given (for two surfaces 
of length dDNA in the x direction) by [11] 



C^DNA 



dh (x) 
dx 



dx. 



(29) 



The crossover from long scale bending-dominated to short scale surface tension-dominated elasticity occurs on length 
scale I which is of the order of I ~ where n is the bending modulus. Our previous measurements [s^l give 

K ~ AOksT and / ~ 50 A, which yields 7^ ~ 25 x 10~^ J/m^. Values of 7p measured for other computer models of 
bilayer membranes (sol . [5l| have been found to be of the same order of magnitude - which happens to be comparable 
to the surface tension of water-oil interfaces. Notice that the superscript "1" rather than "2" is used in Eq. (|29p . 
i.e., we consider this term as linear in the expansion parameters despite of the fact that the elastic energy appears 
as a second order term in (/i/dDNA)- This is because the surface tension 7^ is an order of magnitude larger than the 
energy per unit area cr^,/dDNA/e appearing in (j24p and, therefore, these two terms are comparable to each other. 
The surface tension 7^ is also an an order of magnitude larger than the thermal energy per lipid area fc^T/aiipid, 
appearing in the second order term (j28p which is associated with the mixing entropy. 

With the above expressions for the various terms in the free energy of the system, the equilibrium profile h{x) of the 
surfaces is found by minimizing the leading first order contribution F^ — F^ + F^. The function h(x) is determined 
by the Euler-Lagrange equation 



d^hjx) 
dx'^ 



0, 



(210) 



where the length scale is ^ = e7p/(27r(Tj^). The solution to Eq. (|210p . with boundary conditions h{0) — /i(c?dna) 
is 

1 



h{x) = —X (dDNA - x) 



(211) 



The maximum deformation of the surfaces is obtained at the center of the unit cell {x — dxmKl'^")- For typical values 
of the physical parameters in the system 



Olipid 



70A^ ; o-M ^ 0.8- 



aiipid 



doNA 30A; 7p - 25 X 10"^ jL, 



(212) 



^(c^dna/2) ^ lA<^ c^DNAj which justifies treating {Jijdxmp^ as a small parameter and explains the apparent flatness 
of the surfaces observed in the simulations (e.g., Fig. [T|). It should be noted that our analysis ignores the very small 
overlap that exists between the surface h{x) and the excluded volume of the DNA rods near the edges of the unit cell. 

The charge density fluctuations ba(x) will be determined by minimizing the free energy given by the sum of second 
order terms = F'^ -I- F| -I- -I- F4 , under the total charge conservation constraint (|^^ and with h[x^ = h{x). 
Expressing 6a(x) as a Fourier series 



4- 00 



2-!iin{x/duNA) 



(213) 



we find, after some algebra, that the optimal charge distribution is obtained for the Fourier series 



Cn. — 



CA/^DNA 



ksT/ (aiipidO-A/ [e/aiipid - (Tm]) + c^dna/ {e\n\) 



(214) 



where S in the numerator denotes Kronecker's delta. An approximate, but more useful, form for Sa{x) can be obtained 
by noting that typical values of the physical parameters [see Eq. (|212p ] lead to the first term in the denominator being 
larger than the second term for all values of n. Thus, a reasonable approximation may be obtained by dropping the 
smaller term, which effectively means neglecting the term F| [Eq. ((27| ] in the second order free energy. Without this 
term, the real space form of Sa is given by 



Sa{x) 



CA/aiipidX Z^DNA^S „ ( f 27rx \ 

' B < COS + A 

fllipid L V "DNA / 



X (d 



DNA 



d^ 

"DNA 



(215) 
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where Ib = e^/(e/cBT) 7.lA is the Bjerum length, and A = B~^{it/2) ■ (doNA/O- From this expression we readily 
conclude the following: (1) Only for A — 0, which corresponds to the limit of infinitely rigid surfaces (7p oo), is 
the peak of the charge distribution at the edge of the unit cell. (2) Upon increasing A, the maximum of 5a{x) shifts 
gradually towards the center of the unit cell, and (3) stays in the center, x = doNA/Z, for all values A > 27r^. One 
can easily verify that for the range of parameters in our simulations, < ^ < 27r^, which explains our observation 
of the maximum of 6a{x) somewhere between the edge and the center of the unit cell (see Fig. [T]). We can also use 
Eq. (|215p to estimate the amplitude of the charge fluctuations. For the physical parameters given in Eq. (|212p . we 
have B ^ 0.12. The maximum of 5a{x) is obtained for x ^ O.ldDNA (which should be compared to a; 0.2(iDNA in the 
simulations), where Sa/aM ~ 0.07 (compare to Sa/am 0.1 in the simulations). We consider this semi-quantitative 
agreement with the numerical results as reasonable, given the approximate nature of our analytical model. 
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FIG. 1: Equilibrium configuration of a complex consisting of two bilayer membranes, each with 390 lipids and five DNA strands. 
The lipids are modeled as trimers with hydrophilic (black) and hydrophobic (gray) particles. The DNA (red) are modeled as 
rigid rods with a uniform negative axial charge density. Then complex is isoelectric, i.e., the negative charge of the DNA is 
neutralized by the charge of the cationic lipids with no added salt. Thus, each bilayer in the shown complex includes 150 
monovalent lipids, all of which reside in the inner layers facing the DNA array. 
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FIG. 2: Average DNA spacing, doNA, as a function of the inverse of the fraction of charged hpids l/<j>c- Markers - numerical 
results (uncertainties are smaller than symbols); solid line - fit to Eq. ([S]) with aupid = 69A^. 
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FIG. 3: The average area per lipid, aupid, as a function of the fraction of charged hpids 4>c- 
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FIG. 4: The effective stretching modulus, K^, of the complex as a function of 0c- 



FIG. 5: Equilibrium configuration of a complex with (j>c ~ 0.9 whose membranes develop pores. 




FIG. 6: Local fraction of charged lipids (j>c as a function of x, the position within a unit cell of the complex. Curves, from 
bottom to top, correspond to mean fraction of 0.3, 0.4, 0.5, 0.6, 0.68, 0.77, 0.81. 
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FIG. 7: Local chaxge density of the membranes <tm as a function of x. Curves, from bottom to top, correspond to 0c ~ 0.3, 0.4, 
0.5, 0.6, 0.68, 0.77, 0.81. Dashed horizontal line corresponds to the effective charge density of the DNA ctdna ~ 9.4 •lO-^e/A^ 




FIG. 8: Total area density of the lipids p as a function of x for different values of 0c- 
density of uncharged membranes. Lines are guide to the eyes. 
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= 1/69A-^ is the area 




FIG. 9: Schematic picture of the complex consisting of an array of equally spaced DNA rods with nearest neighbor spacing 
dDNA and two surfaces separated a distance D from the mid-plane of the DNA array. The DNA rods are uniformly charged 
with charge density A < per unit length. Surfaces have a mean charge density ctm > per unit area and local charge density 
am + 5a{x). Their local height above/below the DNA mid-plane is denoted by D — h{x). Lower surface is drawn in the reference 
state, where 5a = 0, and, h = Q. 



